#include "frequency.h"
#include <iostream>
#include <fstream>
#include <cstring>
#include <cmath>

using namespace std;

/*************************************************************************
* Profile -  constructor
*   parameters:
*       -kmer_size: k
*       -r: read size
*       -x: file name of the x axis reference sequence
*       -y: file name of the y axis reference sequence
*   all other class variables are set to zero or NULL; axis are loaded
**************************************************************************/
Profile::Profile(unsigned short int kmer_size, int r, char * x, char * y)
{
    int i;
    k=kmer_size;
    size=pow(4,k);
    sequence=NULL;
    sequence_len=0;
    read_size = r;
    frequencies = new unsigned int [size];
    tableX = new unsigned int [size];
    tableY = new unsigned int [size];
    for(i=0; i<size; i++) {frequencies[i]=0; tableX[i]=0; tableY[i]=0;}             //initialize all values to 0

    load_axis(x,0);                             //makes first file "char * x" the x axis
    load_axis(y,1);                             //makes first file "char * y" the y axis

    circular=false;
}

/*************************************************************************
* add_read
*   parameters:
*       -read_seq:  the read sequence as a character
*   returns true if the read is of the right size and read_seq contains
*   a sequence
**************************************************************************/
bool Profile::add_read(char * read_seq)
{
    if(sequence!=NULL)
        delete sequence;
    sequence_len=read_size;
    sequence=new char [sequence_len+1];
    strcpy(sequence,read_seq);
    return true;
}

/*************************************************************************
* get_sequence -  gets sequence from file and stores in sequence
*   parameters:
*       -file_name: the file with the FASTA format sequence
*   return: true if add successful, false if otherwise
**************************************************************************/
bool Profile::get_sequence(char * file_name)                                                              //WORKS -- TESTED 120501 CP
{
    char c;
    ifstream in;
    int i;
    char header[500];
    sequence_len=0;
    if(sequence!=NULL)
        delete sequence;
    in.open(file_name);
    if (!in.is_open())                                   //file is in FASTA format
    {
        cout << "sequence not opened" << endl;
        return false;
    }
    in.getline(header,500,'\n');
    while(in.peek()!=EOF)
    {
        in>>c;
        sequence_len++;
    }
    in.clear();
    in.close();

    sequence=new char [sequence_len+1];

    in.open(file_name);
    in.getline(header,500,'\n');
    for(i=0; i<sequence_len; i++)
    {
        in>>sequence[i];
        if(sequence[i]>96) sequence[i]-=32;
    }
    sequence[i]='\0';
    in.clear();
    in.close();

    return true;
}

/*************************************************************************
* Profile -  deconstructor
*   delete all dynamically allocated class variables
**************************************************************************/
Profile::~Profile()
{
    delete frequencies;
    delete tableX;
    delete tableY;
    delete sequence;
}

/*************************************************************************
* examine_read
*   parameters:
*       -read_seq: the read sequence
*       -c: the coordinates within the library space returned to user
*   returns true if able to examine the read, else false
**************************************************************************/
bool Profile::examine_read(char * read_seq, coordinates &c)
{
    if(add_read(read_seq))
    {
        if(get_coordinate_read(c)) return true;
        else return false;
    }
    else return false;
}

/*************************************************************************
* get_coordinate_read
*   parameters:
*       -c: coordiante to be returned to user for sequence
*   Makes calls to determine where read falls in k-mer space
*   returns true if the sequence is loaded correctly.
**************************************************************************/
bool Profile::get_coordinate_read(coordinates &c)
{
    if(sequence==NULL || sequence_len<read_size) return false;
    calc_freqs(0,read_size-1);
    c.x=correlation(0);         //X axis
    c.y=correlation(1);         //Y axis
    return true;
}


/*************************************************************************
* get_coordinates_genome -  gets all of the coordinates for the reference
*                           genome sequence
*   parameters:
*       -file_name: the file listing all of the reference genomes
*       -refr_num:  the reference genome # in the file
*       -c:         the array of coordinates to be populated
*   return: number of coordinates
**************************************************************************/
unsigned int Profile::get_coordinates_genome(char * file_name, int refr_num, coordinates *& c)            //WORKS -- TESTED 120501 CP
{
    std::vector<float> x;
    std::vector<float> y;
    int num_coordinates;
    int i;
    char line[1000];
    char top[1000];

    ifstream in;

    in.open(file_name);
    if (in.is_open()==false)
    {
        cout << "something is wrong" << endl;
    }
    for(i=0; i<refr_num; i++)
    {
        in.getline(line,1000,'\n');             //gets reference genome information
        while(in.peek()=='-')                   //it's not a reference genome
        {
            if(in.peek()==EOF)
            {
                cout << "ERROR with format of file of reference genomes. Analysis terminated.\n";
                return 0;
            }
            in.getline(line,1000,'\n');
        }
    }

    //loop through reference genome and all of its alternate genomes
    in.getline(line,1000,'\t');             //gets genus species information
    in.getline(line,1000,'\t');             //gets reference file name
    in.getline(line,1000,'\t');             //gets reference file name
    in.getline(top,1000,'\n');             //gets topology
    if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
        circular=true;
    else
        circular=false;
    get_sequence(line);                     //this is the reference
    num_coordinates=get_coordinates_genome(x,y);

    while(in.peek()=='-')
    {
        in.getline(line,1000,'\t');             //gets genus species information
        in.getline(line,1000,'\t');             //gets reference file name
        in.getline(line,1000,'\t');             //gets alternate g/s file name
        in.getline(top,1000,'\n');             //gets topology
        if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
            circular=true;
        else
            circular=false;
        if(sequence!=NULL)
            delete sequence;
        sequence = NULL;
        get_sequence(line);
        num_coordinates+=get_coordinates_genome(x,y);           //will add to existing vectors
    }
    in.clear();
    in.close();


    c=new coordinates [num_coordinates];

    for(i=0; i<num_coordinates; i++)
    {
        c[i].x=x[i];
        c[i].y=y[i];
    }
    x.clear();
    y.clear();

    return num_coordinates;
}

/*************************************************************************
* get_centroids_genome -  gets all of the coordinates for the reference  //MS: changed from coordinates with probability taken as percentage of genome
*                           genome sequence and outputs centroids
*   parameters:
*       -file_name: the file listing all of the reference genomes
*       -refr_num:  the reference genome # in the file
*       -c:         the array of centroids to be populated
*   return: number of centroids
**************************************************************************/
unsigned int Profile::get_centroids_genome(char * file_name, int refr_num, centroid *& c)
{
    std::vector<float> x;
    std::vector<float> y;
    std::vector<float> prob;
    int num_centroids;
    int i;
    char line[1000];
    char top[1000];

    ifstream in;

    in.open(file_name);
    if (in.is_open()==false)
    {
        cout << "something is wrong" << endl;
    }
    for(i=0; i<refr_num; i++)
    {
        in.getline(line,1000,'\n');             //gets reference genome information
        while(in.peek()=='-')                   //it's not a reference genome
        {
            if(in.peek()==EOF)
            {
                cout << "ERROR with format of file of reference genomes. Analysis terminated.\n";
                return 0;
            }
            in.getline(line,1000,'\n');
        }
    }

    //loop through reference genome and all of its alternate genomes
    in.getline(line,1000,'\t');             //gets genus species information
    in.getline(line,1000,'\t');             //gets reference file name
    in.getline(line,1000,'\t');             //gets reference file name
    in.getline(top,1000,'\n');             //gets topology
    if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
        circular=true;
    else
        circular=false;
    get_sequence(line);                     //this is the reference

    num_centroids=get_centroids_genome(x,y,prob);

    while(in.peek()=='-')
    {
        in.getline(line,1000,'\t');             //gets genus species information
        in.getline(line,1000,'\t');             //gets reference file name
        in.getline(line,1000,'\t');             //gets alternate g/s file name
        in.getline(top,1000,'\n');             //gets topology
        if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
            circular=true;
        else
            circular=false;
        if(sequence!=NULL)
            delete sequence;
        sequence = NULL;
        get_sequence(line);
        num_centroids+=get_centroids_genome(x,y,prob);           //will add to existing vectors
    }
    in.clear();
    in.close();


    c=new centroid [num_centroids];

    for(i=0; i<num_centroids; i++)
    {
        c[i].x=x[i];
        c[i].y=y[i];
        c[i].prob=prob[i];
    }
    x.clear();
    y.clear();

    return num_centroids;
}

/*************************************************************************
* get_centroids_genome -  gets all of the coordinates for the reference  //MS: changed to centroids as above, changed to vectors
*                           genome sequence and outputs centroids
*   parameters:
*       -file_name: the file listing all of the reference genomes
*       -refr_num:  the reference genome # in the file
*       -c:         the array of centroids to be populated
*   return: number of centroids
**************************************************************************/
unsigned int Profile::get_centroids_genome(char * file_name, int refr_num, vector<centroid> & centros)            //WORKS -- TESTED 120501 CP
{
//    std::vector<float> x;
//    std::vector<float> y;
//    std::vector<float> prob;
//    int num_centroids;
    int i;
    char line[1000];
    char top[1000];

    ifstream in;

    in.open(file_name);
    if (in.is_open()==false)
    {
        cout << "something is wrong" << endl;
    }
    for(i=0; i<refr_num; i++)
    {
        in.getline(line,1000,'\n');             //gets reference genome information
        while(in.peek()=='-')                   //it's not a reference genome
        {
            if(in.peek()==EOF)
            {
                cout << "ERROR with format of file of reference genomes. Analysis terminated.\n";
                return 0;
            }
            in.getline(line,1000,'\n');
        }
    }

    //loop through reference genome and all of its alternate genomes
    in.getline(line,1000,'\t');             //gets genus species information
    in.getline(line,1000,'\t');             //gets reference file name
    in.getline(line,1000,'\t');             //gets reference file name
    in.getline(top,1000,'\n');             //gets topology
    if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
        circular=true;
    else
        circular=false;
    get_sequence(line);                     //this is the reference

    get_centroids_genome(centros);

    /**
    eliminated with change to single genome paradigm, eliminates need to search through fof
    **/
//    while(in.peek()=='-')
//    {
//        in.getline(line,1000,'\t');             //gets genus species information
//        in.getline(line,1000,'\t');             //gets reference file name
//        in.getline(line,1000,'\t');             //gets alternate g/s file name
//        in.getline(top,1000,'\n');             //gets topology
//        if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
//            circular=true;
//        else
//            circular=false;
//        if(sequence!=NULL)
//            delete sequence;
//        sequence = NULL;
//        get_sequence(line);
//        num_centroids+=get_centroids_genome(x,y,prob);           //will add to existing vectors
//    }
    in.clear();
    in.close();


//    c=new centroid [centros.size()];

//    for(i=0; i<num_centroids; i++)
//    {
//        c[i].x=x[i];
//        c[i].y=y[i];
//        c[i].prob=prob[i];
//    }
//    x.clear();
//    y.clear();

    return centros.size();
}

/*************************************************************************
* get_centroids_genome -  gets all of the coordinates for the reference     //elimination of reference number since 1 genome per file
*                           genome sequence and outputs centroids
*   parameters:
*       -file_name: the file listing all of the reference genomes
*       -refr_num:  the reference genome # in the file
*       -c:         the array of centroids to be populated
*   return: number of centroids
**************************************************************************/
unsigned int Profile::get_centroids_genome(char * file_name, vector<centroid> & centros)
{
//    std::vector<float> x;
//    std::vector<float> y;
//    std::vector<float> prob;
//    int num_centroids;
//    int i;
//    char line[1000];
//    char top[1000];
//
//    ifstream in;
//
//    in.open(file_name);
//    if (in.is_open()==false)
//    {
//        cout << "something is wrong" << endl;
//    }
////    for(i=0; i<refr_num; i++)
////    {
//        in.getline(line,1000,'\n');             //gets reference genome information
////        while(in.peek()=='-')                   //it's not a reference genome
//        {
//            if(in.peek()==EOF)
//            {
//                cout << "ERROR with format of file of reference genomes. Analysis terminated.\n";
//                return 0;
//            }
//            in.getline(line,1000,'\n');
//        }
//    }

    //loop through reference genome and all of its alternate genomes

//    in.getline(line,1000,'\n');             //gets topology
//    if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
        circular=true;
//    else
//        circular=false;
    get_sequence(file_name);                     //this is the reference

    get_centroids_genome(centros);

//    while(in.peek()=='-')
//    {
//        in.getline(line,1000,'\t');             //gets genus species information
//        in.getline(line,1000,'\t');             //gets reference file name
//        in.getline(line,1000,'\t');             //gets alternate g/s file name
//        in.getline(top,1000,'\n');             //gets topology
//        if( (strcmp(top,"circular")==0) || (strcmp(top,"Circular")==0) || (strcmp(top,"CIRCULAR")==0) )
//            circular=true;
//        else
//            circular=false;
//        if(sequence!=NULL)
//            delete sequence;
//        sequence = NULL;
//        get_sequence(line);
//        num_centroids+=get_centroids_genome(x,y,prob);           //will add to existing vectors
//    }
//    in.clear();
//    in.close();


//    c=new centroid [centros.size()];

//    for(i=0; i<num_centroids; i++)
//    {
//        c[i].x=x[i];
//        c[i].y=y[i];
//        c[i].prob=prob[i];
//    }
//    x.clear();
//    y.clear();

    return centros.size();
}

/*************************************************************************
* get_coordinates_genome -  gets all of the coordinates for the sequence
*                           PRIVATE -- called by public get_coords funct
*   parameters:
*       -x: vector of correlations for axis X
*       -y: vector of correlations for axis y
*   return: number of coordinates
**************************************************************************/
unsigned int Profile::get_coordinates_genome(std::vector<float> & x, std::vector<float> & y)            //WORKS -- TESTED 120501 CP
{
    int i;
    int counter=0;
    if(!circular)
    {
        for(i=0; i<(sequence_len-read_size+1); i++)
        {
            calc_freqs(i,i+read_size-1);
            x.push_back(correlation(0));         //X axis
            y.push_back(correlation(1));         //Y axis
            counter++;
        }
    }
    else                                         //need to calculate frequencies over the origin
    {
        for(i=0; i<sequence_len; i++)
        {
            calc_freqs(i,i+read_size-1);
            x.push_back(correlation(0));         //X axis
            y.push_back(correlation(1));         //Y axis
            counter++;
        }
    }
    return counter;
}

/*************************************************************************
* get_centroids_genome -  gets all of the centroids for the sequence            //MS: changed from coordinates probability is generated based on % genome reads make up at a point
*                           PRIVATE -- called by public get_coords funct
*   parameters:
*       -x: vector of correlations for axis X
*       -y: vector of correlations for axis y
*   return: number of coordinates
**************************************************************************/
unsigned int Profile::get_centroids_genome(std::vector<float> & x, std::vector<float> & y, std::vector<float> & prob)
{
    int i;
    int counter=0;
    float denom = sequence_len-read_size+1;
    float tempProb = 1/denom;
    if(!circular)
    {
        for(i=0; i<(sequence_len-read_size+1); i++)
        {
            calc_freqs(i,i+read_size-1);
            x.push_back(correlation(0));         //X axis
            y.push_back(correlation(1));         //Y axis
            prob.push_back(tempProb);
            counter++;
        }
    }
    else                                         //need to calculate frequencies over the origin
    {
        for(i=0; i<sequence_len; i++)
        {
            calc_freqs(i,i+read_size-1);
            x.push_back(correlation(0));         //X axis
            y.push_back(correlation(1));         //Y axis
            prob.push_back(tempProb);
            counter++;
        }
    }
    return counter;
}

/*************************************************************************
* get_centroids_genome -  gets all of the centroids for the sequence            //MS: simplified version from above, merge x, y, prob into cent struct
*                           PRIVATE -- called by public get_coords funct
*   parameters:
*       -x: vector of correlations for axis X
*       -y: vector of correlations for axis y
*   return: number of coordinates
**************************************************************************/
unsigned int Profile::get_centroids_genome(vector<centroid> & c)
{
    int i;
    int counter=0;
    float denom = sequence_len-read_size+1;
    float tempProb = 1/denom;
    if(!circular)
    {
        for(i=0; i<(sequence_len-read_size+1); i++)
        {
            centroid tempC;
            calc_freqs(i,i+read_size-1);
            tempC.x = correlation(0);         //X axis
            tempC.y = correlation(1);         //Y axis
            tempC.prob = tempProb;
            c.push_back(tempC);
            counter++;
        }
    }
    else                                         //need to calculate frequencies over the origin
    {
        for(i=0; i<sequence_len; i++)
        {
            centroid tempC;
            calc_freqs(i,i+read_size-1);
            tempC.x = correlation(0);         //X axis
            tempC.y = correlation(1);         //Y axis
            tempC.prob = tempProb;
            c.push_back(tempC);
            counter++;
        }
    }
    return counter;
}

/*************************************************************************
* calc_freqs -  calculates frequencies of k-mers in sequence
*               updates class variable frequencies
*   parameters:
*       -start: start position in sequence for k-mer calculation
*       -end: end position in sequence for k-mer calculation
*   if going over the origin, end>sequence_len && circular=true
*   return: true if add successful, false if otherwise
**************************************************************************/
bool Profile::calc_freqs(unsigned int start, unsigned int end)                                                              //WORKS -- TESTED 120501 CP
{
    int i, j;
    unsigned int value;
    bool oddchar;

    if(start<0)         //if start is out of bounds
        return false;
    if((end>sequence_len) && (circular==false))
        return false;

    int o_end=0;
    if(end>sequence_len) {o_end=end; end=sequence_len-1;}

    for(i=0; i<size; i++) frequencies[i]=0;
   // cout << "1 k " << k << " end " << end << "end-k+1 " << end-k+1 << endl;
    for(i=start; i<=(end-k+1); i++)
    {
        //cout << "2 k " << k << " end " << end << "end-k+1 " << end-k+1 << endl;
        oddchar=false;
        value=0;
        //cout << "sequence: " << sequence << endl;
        for(j=0; j<k; j++)
        {
            //cout << " i " << i << " j " << j << " i + j " << i+j << " seq length: " << sequence_len << " k " << k << endl;
            switch(sequence[i+j])
            {
            case 'A':
                value=value << 2;
                break;
            case 'T':
                value=value << 2;
                value += (int) 1;
                break;
            case 'C':
                value=value << 2;
                value += (int) 2;
                break;
            case 'G':
                value=value << 2;
                value += (int) 3;
                break;
            default:
                oddchar=true;
                j=k;
                break;
            }
        }
        if(!oddchar)
            frequencies[value]++;
    }
    //reverse complement of sequence
    for(i=start; i<=(end-k+1); i++)
    {
        oddchar=false;
        value=0;
        for(j=k-1; j>=0; j--)
        {
            switch(sequence[i+j])
            {
            case 'T':
                value=value << 2;
                break;
            case 'A':
                value=value << 2;
                value += (int) 1;
                break;
            case 'G':
                value=value << 2;
                value += (int) 2;
                break;
            case 'C':
                value=value << 2;
                value += (int) 3;
                break;
            default:
                oddchar=true;
                j=0;
                break;
            }
        }
        if(!oddchar)
            frequencies[value]++;
    }


    if(circular && (o_end!=0))                                                          //account for going over the origin -- circular sequence
    {
        int size_temp=(k-1)+(o_end%sequence_len);
        char * temp_sequence=new char [size_temp];
        j=0;
        for(i=sequence_len-k+1; i<sequence_len; i++)
        {
            temp_sequence[j]=sequence[i];
            j++;
        }
        for(i=0; i<(size_temp-(k-1)); i++)
        {
            temp_sequence[j]=sequence[i];
            j++;
        }

        for(i=0; i<(o_end%sequence_len); i++)                     //k-mers over origin
        {
            oddchar=false;
            value=0;
            for(j=0; j<k; j++)
            {
                switch(temp_sequence[i+j])
                {
                case 'A':
                    value=value << 2;
                    break;
                case 'T':
                    value=value << 2;
                    value += (int) 1;
                    break;
                case 'C':
                    value=value << 2;
                    value += (int) 2;
                    break;
                case 'G':
                    value=value << 2;
                    value += (int) 3;
                    break;
                default:
                    oddchar=true;
                    j=k;
                    break;
                }
            }
            if(!oddchar)
                frequencies[value]++;
        }

        //reverse complement of sequence
        for(i=0; i<(o_end%sequence_len); i++)
        {
            oddchar=false;
            value=0;
            for(j=k-1; j>=0; j--)
            {
                switch(temp_sequence[i+j])
                {
                case 'T':
                    value=value << 2;
                    break;
                case 'A':
                    value=value << 2;
                    value += (int) 1;
                    break;
                case 'G':
                    value=value << 2;
                    value += (int) 2;
                    break;
                case 'C':
                    value=value << 2;
                    value += (int) 3;
                    break;
                default:
                    oddchar=true;
                    j=0;
                    break;
                }
            }
            if(!oddchar)
                frequencies[value]++;
        }
        delete temp_sequence;
    }


    return true;
}

/*************************************************************************
* correlation -  calculates the correlation between frequency & axis
*   parameters:
*       -axis: 0 for frequencies loaded as X axis, 1 for Y axis
*   return: float correlation value scaled between 0 and 1
**************************************************************************/
float Profile::correlation(int axis)                                                              //WORKS -- TESTED 120501 CP
{
    float xy=0;
    float x=0;
    float y=0;
    float xx=0;
    float yy=0;

    for (int i=0; i<size; i++)
    {
        if(axis==0)
        {
            xy+=(frequencies[i]*tableX[i]);
            y+=tableX[i];
            yy+=pow(tableX[i],2);
        }
        else
        {
            xy+=(frequencies[i]*tableY[i]);
            y+=tableY[i];
            yy+=pow(tableY[i],2);
        }
        x+=frequencies[i];
        xx+=pow(frequencies[i],2);
    }

    float num= xy-((x*y)/size);
    float Square=((xx-((x)*(x))/size)*(yy-((y)*(y))/size));
    float sqr= sqrt(Square);
    float R= num/sqr;

    //scale r between 0 and 1 (to avoid negative values)
    if(R==-1.0) return 0;
    else
        R=(R+1.0)/2.0;

    return R;
}

/*************************************************************************
* load_axis -  calculates the correlation between frequency & axis
*   parameters:
*       -file_name: name of the reference GC file
*       -axis: 0 for frequencies loaded as X axis, 1 for Y axis
*   return: true if loaded, else false
**************************************************************************/
bool Profile::load_axis(char * file_name, int axis)                                                              //WORKS -- TESTED 120501 CP
{
    if(axis>1) return false;

    ifstream in;
    in.open(file_name);
    int v,value;
    char c;
    int i;

    in>>c;
    in>>v;

    while(v!=k)
    {
        for(i=0; i<pow(4,v); i++)
            in>> value;
        in>>c;
        in>>v;
    }
    if(axis==0)
    {
        for(i=0; i<size; i++)
            in>> tableX[i];
    }
    else
    {
         for(i=0; i<size; i++)
            in>> tableY[i];
    }
    in.clear();
    in.close();
    return true;
}
